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ABSTRACT 

MicroRNA (miRNA) target hubs are genes that can 
be simultaneously targeted by a comparatively large 
number of miRNAs, a class of non-coding RNAs that 
mediate post-transcriptional gene repression. 
Although the details of target hub regulation 
remain poorly understood, recent experiments 
suggest that pairs of miRNAs can cooperate if 
their binding sites reside in close proximity. To 
test this and other hypotheses, we established a 
novel approach to investigate mechanisms of col- 
lective miRNA repression. The approach presented 
here combines miRNA target prediction and tran- 
scription factor prediction with data from the litera- 
ture and databases to generate a regulatory map for 
a chosen target hub. We then show how a kinetic 
model can be derived from the regulatory map. To 
validate our approach, we present a case study for 
p21, one of the first experimentally proved miRNA 
target hubs. Our analysis indicates that distinctive 
expression patterns for miRNAs, some of which 
interact cooperatively, fine-tune the features of tran- 
sient and long-term regulation of target genes. With 
respect to p21 , our model successfully predicts its 
protein levels for nine different cellular functions. In 
addition, we find that high abundance of miRNAs, in 
combination with cooperativity, can enhance noise 
buffering for the transcription of target hubs. 



INTRODUCTION 

The regulation of basic cellular functions is controlled by 
complex networks involving interacting genes, proteins 
and small molecules (1). In the past decade, an additional 
level of regulation was discovered after the identification 
of a class of short non-coding RNAs called microRNAs 
(miRNAs). miRNAs are ~22nt in length and their main 
function is to regulate the activity and stability of specific 
mRNAs at the post-transcriptional level (1,2). Individual 
miRNAs can regulate a large number of mRNA tar- 
gets and target genes can be regulated by multiple 
miRNAs (3,4). 

The term target hub was introduced by Borneman et al. 
(5), who constructed a putative transcriptional control 
network for yeast cell differentiation. Genes, combining 
several transcription factors (TFs) in their promoter 
regions, were referred to as target hubs. Shortly thereafter, 
Shalgi et al. (3) defined the notion of miRNA target hub 
for genes that are regulated by at least 15 different 
miRNAs (Figure 1). 

At least two groups (6,7) experimentally proved that 
pairs of miRNAs, with binding sites that reside in close 
proximity or partially overlap, show cooperative behav- 
iour, adding to the complexity of the miRNA target regu- 
lation. Saetrom et al. (7) proposed a range of 13-35 nt for 
the distance of the binding start sites for which the 
phenomenon appears. According to these experimental 
evidences, we assumed this kind of interaction as plausible 
and applied it in a kinetic model of p21 regulation by 
multiple miRNAs to investigate the dynamic and 
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Figure 1. Illustration of miRNA target hub regulatory network. Once an miRNA target hub gene is transcribed, its mRNAs are 
post-transcriptionally regulated by many miRNAs whose expression can be activated or inhibited by numerous TFs. In animals, the complementarity 
between an miRNA and the target hub mRNA usually results in reduced protein expression through translational repression or mRNA 
deadenylation followed by degradation. Target hubs are highly interconnected via protein-protein interactions, so that its regulatory network 
may contain interlocked network motifs, such as feedback and feedforward loops. 



regulatory consequences of this feature. A well-known 
miRNA target hub is the cell cycle regulator and 
tumour suppressor p21 (a.k.a. CDKN1A or Cipl/Wafl) 

(8) . Specifically, p21 is involved in the Gi phase cell cycle 
arrest in response to stress signals provoked by DNA 
damage. The expression of p21 is dependent on environ- 
mental conditions and is transcriptionally regulated 
through p53-dependent and -independent mechanisms 

(9) . Recently, numerous miRNAs were identified to 
regulate the expression of p21. From 266 predicted 
miRNA regulators of p21, Wu et al. (10) validated the 
ability of a subset of 28 miRNAs to target p21 using a 
high-throughput luciferase reporter screen. Thus, p21, to 
our knowledge, is one of the first experimentally verified 
miRNA target hubs. Moreover, the regulation induced by 
miRNAs can provoke phenotypic changes. Wu et al. (10) 
showed that overexpression or knockout of miR-372, 
miR-519-3p and miR-520a-3p in JAR cells leads to 
either increased cell proliferation or accumulation of 
cells in Gl phase. 

In this work, we investigated the structure and 
dynamics of the regulatory network of the miRNA 
target hub p21 and several of its repressor miRNAs by 
using a quantitative kinetic model. By integrating infor- 
mation and data from the literature and databases, we 
constructed a regulatory map of p21 including its target- 
ing miRNAs, TFs and interacting proteins. Using this 
map, we found potential 'housekeeping' miRNAs that 
are promoted by a few TFs and are associated with 
many distinct cellular functions. Based on this map, we 
derived a kinetic model to test hypotheses about mechan- 
isms of collective miRNA repression. We computationally 



detected several couples of miRNAs that may cooperate in 
the repression of p21 and experimentally validated the 
cooperativity between two of them, miR-572 and 
miR-93. We further used these two miRNAs as an 
example to investigate three hypothetical target regulation 
mechanisms by multiple miRNAs. We found that one 
miRNA is sufficient to completely repress a target when 
highly expressed. However, the same effect can be 
achieved by two cooperating miRNAs, both of which 
are only expressed on a medium level. The model was 
further used to predict miRNA and p21 expressions for 
nine cellular functions, and it predicted high levels of p21 
expression during DNA damage, DNA repair, senescence 
and migration, whereas low expressions of p21 for cell 
proliferation and apoptosis. 

MATERIALS AND METHODS 

General overview of the methodology 

The aim of our analysis is to unravel the complex mech- 
anisms by which gene and signalling networks are 
regulated by multiple miRNAs. We iteratively integrate 
data from the biomedical literature, high-throughput ex- 
periments and biological databases into a kinetic model of 
miRNA target hub regulation. The kinetic model is then 
used to formulate and test hypotheses about mechanisms 
of target regulation and cellular function-related vari- 
ability. In short, the approach includes three blocks 
(Figure 2): (i) data retrieval, (ii) construction and 
analysis of the regulatory map and (iii) kinetic modelling 
and simulation. These blocks are discussed in detail below. 
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Figure 2. Schematic representation of the approach. Biomedical knowledge and quantitative data are integrated using mathematical modelling to 
describe the regulatory role of miRNAs in signalling pathways and gene regulatory networks. It considers data retrieval, iterative cycles of model 
construction and calibration and computational simulations. The mathematical model obtained is suitable to test hypotheses and simulate complex 
biological scenarios associated with cellular function-related variability. 



Data retrieval 

Information about protein interactions was extracted 
from the Human Protein Reference Database (release 
9.0) (11) and the database STRING (release 9.0) (12). A 
list of experimentally verified TFs for the considered target 
hub was generated from the literature. This list was com- 
plemented by putative TFs that are associated with 
conserved TF-binding sites (human, mouse and rat) 
residing in the 5-kb upstream region of the target gene. 
This information was extracted from the table of TFs with 
conserved binding sites (hgl8) in the USCS genome 
browser (13). Information about miRNA-target inter- 
actions can be extracted either from databases of validated 
interactions, e.g. miRecords (14), Tarbase (15) and 
miRTarBase (16), or from predicted ones, which can be 
found in databases like miRWalk (17) and miRGen 2.0 
(18). In the case of p21, we used the information from Wu 
et al. (10) where a list of predicted p21 -regulating miRNAs 
was subjected to experimental validation. There exist 
other miRNAs than those included in our network that 
can regulate p21 expression (10,19). However, in order to 
enable the construction of a kinetic model, we considered 
only those miRNAs for which there are consistent quan- 
titative data describing their repression abilities on p21. 



A list of TFs controlling the expression of the miRNAs 
was constructed using information of experimentally 
proved TFs of miRNAs contained in TransmiR (release 
1.0) (20). In addition, we generated a list of putative TFs 
of miRNAs with binding sites in the 10-kb upstream 
region of the miRNA genes from the databases PuTmiR 
(release 1.0) and MIR@NT@N (version 1.2.1) and from 
the table of TFs with conserved binding sites 
(tfbsConsFactor) in USCS hgl8 (13,21,22). We selected 
this region following the method in Shalgi et al. (3), 
where the 10-kb upstream region is indicated as 'putative 
regulatory region of miRNAs'. For further investigation, 
we considered only the miRNA TFs that are experimen- 
tally verified or predicted in two of three resources. 

Construction and analysis of the regulatory map 

We first converted the information described in the 
previous section into a regulatory network map in 
CellDesigner (23). For assessing the reliability of the 
network structure, we computed a confidence score for 
each interaction included in the p21 regulatory map, 
similar to the procedure used in several interaction data- 
bases like STRING, iRefWeb or IntAct. The computed 
scores ranged between 0 (totally uncertain interaction) and 
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1 (most reliable), and they were established by integrating 
weighted scores for publications reporting an interaction, 
experimental method(s) used, interaction type and compu- 
tational predictions. The scoring system was inspired by 
the one used in IntAct (24). A detailed description of how 
the confidence scores were calculated is provided in 
Supplementary Materials Section 1. All the interactions 
(p21-protein, TF-p21, TF-miRNA and miRNA-p21) 
displayed in the map and their corresponding scores can 
be found in Supplementary Excel file. 

In our map, the network motifs including a TF that 
regulates p21 and a p21 -targeting miRNA were considered 
as feedforward loops. Those in which p21 is repressed by 
a TF and a targeting miRNA were classified as coherent 
feedforward loops. Those in which p21 is activated by a 
TF and repressed by an miRNA were classified as inco- 
herent feedforward loops. Meanwhile, the network motifs 
including a TF that promotes a targeting miRNA and 
interacts with p21 were classified as feedback loops. 
Based on Gene Ontology (GO) terms, we derived 
tailored TF lists of p21 and its targeting miRNAs for dif- 
ferent cellular functions (cell proliferation, apoptosis, 
immune response, inflammation response, cell cycle 
control, DNA damage response, cell senescence, DNA 
repair and cell motility and migration; Supplementary 
Excel file). 

Kinetic modelling 

Based on the regulatory map, a kinetic model was con- 
structed using ordinary differential equations (ODEs) 
(25). Precisely, the kinetic model accounts for the evolu- 
tion in time of the mRNA (mp21) and protein (p21) ex- 
pression levels of the miRNA target hub p21, the 
p21-targeting miRNAs considered (miR ir i= 1 ... 15) and 
the complexes formed by targeted mRNA and miRNA, 
[mp21\miR^. In all, the model is constituted by 32 time- 
dependent variables and 64 parameters: 

d 

— mp2\ = k synjm p2\ ■ F act (TF mp 2\) 

- mp2\ ■ (kdeg_mp2i + ^ kqssjnm ■ rniR,) 

i 

d 

— [mp2\\miRj] = k ass ^„ iRi ■ mpl\ • miRi 

- kdegjcpxt ■ [mp2l\miRi], i = 1...15 

— miRj — ksmjniR ■ F act (TF mi x) 
dt 

- m iRj ■ (k^g^ix. - kassjniK ■mp2\-miR i ,i=\...\5 

d 

— pll = k svn _ p2 i ■ mp2\ - kdeg^pzx -p2l 

For the p21 mRNA (mp2T), processes considered in the 
model are basal synthesis (k syn mp 2i) mediated by its TF 
(F acr (TF mp2 i)), basal degradation (k deg _ mp21 ) and associ- 
ation with an miRNA (k ass mi Ri). For each p21 -targeting 
miRNA (miRj), processes considered in the model are 
basal synthesis (k syn miRi ) mediated by its TF 
(F act {TF miR ff), basal degradation (k deg _ miB ]) and associ- 
ation with the p21 mRNA target (k ass mi Rj). For each 



complex ([mp21\miRi\), processes considered in the 
model are association of miRj and p21 mRNA to form a 
complex {k ass _, niR i) and degradation (k deg _ CpX ^). For p21 
protein (p21), processes considered in the model are 
mRNA-mediated synthesis of protein (k syn _p rot ) and deg- 
radation (k<te g _p2i)- Additional algebraic equations 
accounting for the total measurable amounts of p21 
mRNA (iiip21 total) and each p21-targeting miRNA 
(miRj total) were also included 

mp2\roTAL — mp2\+ [mp2\ \miRj] 

i 

tniRtjoTAL = miRj+ [mp2\\miRj\ 

i 

The model variables were normalized to the 
non-repressed levels of p21 mRNA and protein (mp21 
and p21 levels were assumed to be 1 when no miRNA 
repression occurs). Values of the model parameters 
characterizing the reaction rates were assigned either by 
using fixed values taken from relevant publications 
(e.g. protein, mRNA and miRNA half-lives) or by 
applying data fitting techniques (26,27). In the second 
case, we processed and used quantitative data, describing 
the repression exerted by the different miRNAs on p21 
mRNA and protein levels (10), to estimate the values of 
the parameters k ass _ m a{i and k deg _ Cp xi- These parameters 
were numerically calculated in biologically meaningful 
intervals using an iterative estimation method, which 
combines global and local optimization algorithms. The 
parameter estimation was carried out using SBtoolbox2 
(28). Computational simulations were performed using 
MATLAB (Mathworks, MA, USA). The construction 
and calibration of the kinetic model is described in detail 
in Supplementary Materials Section 2. 

Experiment methods 

For validating the predictive ability of our kinetic 
model, we performed additional experiments. Precisely, 
SK-Mel-147 melanoma cells were transfected with 
mature miRNA mimics of two miRNAs targeting p21 
(miR-572 and miR-93), either individually at a concentra- 
tion of 100 nM or in combination at 50 nM concentration. 
Cells were pulse-treated with 250 nM of the genotoxic sub- 
stance doxorubicin and protein lysates were collected at 
0, 2, 4, 6, 8, and 24 h after treatment. Expression levels of 
p21 were measured using immunoblotting as explained in 
Supplementary Materials Section 4. 

RESULTS 

The regulatory map of p21 

According to the approach described in the previous 
section, biomedical knowledge about p21 regulation 
retrieved from publications and databases was integrated 
into a regulatory map (Figure 3). The map constructed 
in CellDesigner provides an up-to-date summary of infor- 
mation concerning transcriptional and post- 
transcriptional regulation of p21 and its regulating 
miRNAs. Using CellDesigner, the regulatory map is 
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Figure 3. The regulatory map of p21. Visualization of the annotated SBML file integrating information about the regulation of miRNA target hub 
p21. miRNAs are classified according to their mechanisms of p21 regulation (deadenylation denoted by blue box and translation repression denoted 
by green box). TFs are divided into three groups: TFs regulating only p21 expression (red box), TFs regulating only miRNAs (yellow box) and those 
that regulate the expression of p21 and some of the miRNAs (light blue box). The latter accounts for putative feedforward loops. p21 protein 
interaction partners are clustered into groups (grey boxes). The small boxes, highlighted in green, account for putative feedback loops formed by p21, 
given miRNAs and the interacting partner (e.g. p21->AKTl^miR-132Hp21). Highlighted in blue are TFs that down-regulate expression of the 
considered gene targets. The purple box accounts for cellular functions assigned to the considered TFs. A high resolution version of the regulatory 
map is available in Supplementary Data. 



compliant with the Systems Biology Graphical Notation 
(29). The annotated SBML file of the regulatory map is 
freely available at www.sbi.uni-rostock.de/resources/ 
software/target-hub. Furthermore, the map was comple- 
mented with an annotated Excel file containing accession 



numbers; associated GOs; PubMed identifiers and other 
details of the protein interaction partners, miRNAs and 
TFs. In addition, this file contains the confidence score for 
each interaction in the map. These confidence measures 
were established integrating weighted scores for 



Nucleic Acids Research, 2012, Vol. 40, No. 18 8823 



publications reporting an interaction, experimental 
method(s) used, interaction type and computational pre- 
dictions (see Supplementary Excel file for details). 

Overall, the map suggests a complex regulation of the 
miRNA target hub p21 with multiple transcriptional and 
post-transcriptional regulatory motifs. According to our 
analysis, certain TFs, including EGR1 and SP1, promote 
the expression of groups of miRNAs regulating p21. 
This supports the hypothesis that miRNA cooperativity 
plays a role in the repression of p21. GO analysis shows 
that some of the cellular functions, in particular cell prolif- 
eration, apoptosis and immune response, are prominent 
in the regulation of the p21 -targeting miRNAs by TFs 
(Supplementary Materials Figure S4A). For example, 
miRNA-345 and miR-93 are predicted to be promoted 
by three and seven TFs, respectively. However, the 
coverage of cellular functions of miR-345 (six functions) 
is more extensive than miR-93 (three functions). This 
finding implies that miRNAs like miR-345, which are 
promoted by a few TFs but associated with many cellular 
functions, can be considered as 'housekeeping' miRNA. 
On the other hand, miRNAs like miR-93, which are 
promoted by many TFs but associated with a few cellular 
functions, can be considered cellular function-specific 
miRNAs (Supplementary Materials Figure S4B). 

Kinetic modelling of target hub regulation 
by multiple miRNAs 

We derived a kinetic model with a structure suitable to 
investigate the hypothesis of enhanced miRNA repressive 
ability associated with miRNA-binding sites in close prox- 
imity on the 3'-untranslated region (UTR) of the target 
mRNA (6,7). Enhanced repression in this context origin- 
ates from a cooperative interaction between the two 
miRNAs whose binding sites are close. To test the conse- 
quences of this hypothesis, we first predicted the binding 
sites for p21 -targeting miRNAs on its 3'-UTR and then we 
derived a kinetic model to analyse mechanisms for pairs of 
miRNAs repressing the same target hub. Figure 4 shows 
the matrix of potentially cooperating and non-cooperating 
miRNA pairs targeting p21 for post-transcriptional regu- 
lation (see Supplementary Materials Section 2.1 for ex- 
planation and Supplementary Excel file for binding 
patterns). Our results indicate that several couples of 
miRNAs targeting p21 display the conditions as described 
by the previous studies (6,7) and thus can cause enhanced 
target repression. The only exception is miR-639, whose 
binding site neither overlaps with nor is in close proximity 
to the other miRNA-binding sites. 

Using the approach previously described and based on 
the regulatory map, we derived the kinetic model using 
ODEs, accounting for the evolution in time of the 
mRNA and protein expression levels of the target hub 
p21, the p21 -targeting miRNAs and the complexes 
integrated by the mRNA and the miRNAs. To substanti- 
ate the cooperative effect associated with the proximity of 
miRNA-binding sites, we defined a group of new variables 
([mp21\miRj\miRj\), accounting for ternary complexes 
composed of p21 mRNA and two potential cooperating 
miRNAs (miRi and miR^). For these new variables, we 



considered two processes: (i) the association of p21 
mRNA with miR, and miRj into the complex (k c/axx m mj) 
and (ii) the degradation of the complex (k deg _ CpX i^. Their 
dynamics are represented by the following differential 
equations, which complement the already established set 
of equations for the particular case of miRNA 
cooperativity: 

d 

— mp2\ =k syn _ mp 2\ ■ F ac ,(TF mp 2\) 

- mpll ■ ik deganpll + KsMi ■ miR i 

i 

+ Yl k das S -miRij ■ miR i ' miR j) 
•J 

d , , 

— [mp2l\miRj\miRj] = k dass miRjJ ■ mpll ■ miRi ■ miRj 

~ k deg_c,,Xij ■ [mp2\ \miRi\miRj] 

Model validation 

For validating the structure of the kinetic model, we 
carried out overexpression experiments for miR-572 and 
miR-93, two of the p21 -targeting miRNAs that may co- 
operate to regulate p21 according to our analysis. Cells 
were transfected with those miRNAs either individually 
(100 nM) or in combination (50nM each) and treated 
with doxorubicin, a genotoxic stress-inducing agent. 
Expression levels of p21 were measured by immunoblot- 
ting at different times after doxorubicin treatment. Thus, 
we obtained the data characterizing the dynamics of p21 
expression after genotoxic stress in four scenarios: (i) en- 
dogenous miRNA expression, (ii) only miR-572 
overexpressed (100 nM), (hi) only miR-93 overexpressed 
(100 nM) and (iv) both miRNAs partially overexpressed 
(50 nM each). We further simulated the model by 
configuring it to the designed experimental scenarios and 
compared the model predictions with the experimental 
data (see Supplementary Materials Section 3.1 for 
details). As shown in Figure 5, the model predictions are 
in agreement with the experimental observations for p21 
response after genotoxic stress. Overall, overexpression of 
miR-572 or miR-93 is able to reduce the p21 up-regulation 
after genotoxic stress. In the scenarios 2 and 3 the total 
amount of the transfected miRNA is equal. However, the 
differences in p21 response are due to the different repres- 
sion efficiencies of the two miRNAs. Furthermore, the 
data in scenario 4 indicate that the combined partial 
overexpression of both miRNAs, with binding sites in 
close proximity and therefore potentially cooperating, 
induces stronger down-regulation of p21 response. These 
results validate the hypothesis of cooperativity for this 
couple of miRNAs and also the ability of our method to 
detect and characterize miRNA cooperativity. 

To further test the predictive ability of the model, we 
used it to predict miRNA-regulated p21 expression 
patterns in different biological contexts. First, we ex- 
tracted and normalized the tissue-specific miRNA expres- 
sion levels from the database miRNAmap (release 2.0) 
(30). Next, for each miRNA expression profile obtained, 
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Figure 4. (A) Matrix of cooperating and non-cooperating p21 targeting miRNAs. The intersections of pairs of miRNAs denote their potential 
for cooperation. Grey cells indicate a non-interacting pair and blue cells denote potentially cooperating pairs based on binding site proximity range 
(13-35 nt) defined by Saetrom et al. (7). Red cells, however, indicate that cooperation cannot be established due to binding sites with extensive 
overlap or miRNA pairs that share the same binding site. The figures inside specify the fraction of binding sites interacting with the respective 
partner (the figures in brackets represent the miRNA on the x-axis). (B) Binding sites of regulatory miRNAs in the p21 3'-UTR. 
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Figure 5. Regulation of p21 expression by miR-572 and miR-93. The experimental data of p21 protein expression in response to genotoxic 
stress (Exp) are compared with the model predictions (Sim) in four biological scenarios: (i) both miRNAs are normally expressed (NTC), 
(ii) miR-572 is overexpressed (miR-572), (hi) miR-93 is overexpressed (miR-93) and (iv) both miRNAs are overexpressed (miR-572 + 93). n.u.: 
normalized unit. 
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we computed the steady-state levels of p21 by changing 
accordingly the initial concentrations of the corresponding 
miRNAs in our model. Then, we retrieved qualitative in- 
formation about p21 expression levels in the same tissues 
from the database ArrayExpress (version as of January 
2012) (31). For making the simulation results and experi- 
mental data comparable, the tissue-specific p21 expression 
levels were discretized in two categories (low and high). 
The process is explained in detail in Supplementary 
Materials Section 3.3. Assuming that epigenetic regulation 
exerted by miRNAs plays an important role for tissue- 
specific protein expressions, our model is able to correctly 
predict the relative p21 expression levels in 9 of 12 tissues 
(75%; Supplementary Materials Figure S6). 

Computational analysis of miRNA cooperativity 

Using the calibrated and validated model, we hypothe- 
sized three mechanisms for the target regulation con- 
ducted by pairs of miRNAs: (i) independent target 
regulation, in this case the target site proximity does not 
produce an effective enhanced repression (represented by 
setting k dass mjRlj = 0); (ii) interdependent target regula- 
tion, in which the combined binding of both miRNAs is 
required to repress the target (k uss _ miRi = 0, k ass _ miR j = 0 
and k class miRij > 0) and (iii) synergistic target regulation, 
which is a combination of the previous two target 
regulation mechanisms (k ass miRi > 0, k ass miRj > 0, 
kdass miRij > 0)- We investigated the consequences of the 
three mechanisms for target repression by using the 
miRNA pair miR-93 {miR s ) and miR-572 (miR 14 ) as an 
example (see Supplementary Materials Section 3.1 for 
details). Figure 6 shows the course of the p21 protein ex- 
pression levels for the three repression mechanisms upon 
different initial concentrations of both miRNAs, 
miR 8 (0),miR[ 4 (0) e [10 _1 , 10 2 ]. To make our result com- 
parable with the experimental results generated by Wu 
et al. (10), we computed the expression level of p21 
protein at 48 h. 

In case of independent target regulation (Figure 6A), 
sufficient up-regulation of one of the two miRNAs is 
able to induce the complete repression of p21 protein syn- 
thesis. In combination with the up-regulation of a second 
miRNA, a reduced overall level of the miRNAs is required 
to silence p21. In contrast, effective p21 silencing is 
possible only when both miRNAs are sufficiently 
up-regulated. In other words, the independent overexpres- 
sion of a single miRNA is not able to induce p21 silencing, 
rather results in a low miRNA-specific gene down- 
regulation (Figure 6B). This mechanism may explain the 
relatively poor repression ability of some miRNAs when 
they are experimentally overexpressed alone (4). Finally, 
p21 silencing can be achieved through synergistic target 
regulation in two ways (Figure 6C): (i) independent 
overexpression of any of the two miRNAs, which 
requires high levels of miRNA expression (at least 10 2 - 
fold upregulation) and (ii) combined modest upregulation 
of the two cooperative miRNAs, which reduces the 
amount of miRNA required for gene silencing by at 
least one order of magnitude. Moreover, through this 
mechanism, the system becomes more sensitive to 



miRNA expression changes both for individual and 
combined miRNA down-/up-regulation. 

Based on the previous results (Figure 5), we can say that 
synergistic target regulation is the most likely mechanism 
by which p21 is repressed by the two cooperating miRNAs, 
namely, miR-572 and miR-93. Therefore, we further 
investigated this mechanism for different strengths of the 
miRNA cooperativity (named as K) and simulated four 
scenarios with different initial miRNA concentrations 
([miRg(0), miR 14 (0)]): (i) both miRNAs are normally ex- 
pressed ([1,1], blue triangles), (ii) miR-93 is normally ex- 
pressed, miR-572 is overexpressed ([1,10], red diamonds), 
(iii) miR-93 is overexpressed, miR-572 is normally ex- 
pressed ([10,1], green circles) and (iv) both miRNAs are 
overexpressed ([10,10], purple triangles). For each 
scenario, we defined an interval accounting for the 
strength of the miRNA cooperativity (K = k duss miR s,i4 x 
[10~ 5 10 5 ]). We simulated the model, assigning values of K 
in steps of 10 0 ' 25 ; for every round of simulation, we 
computed the expression levels of p21 protein at 48 h 
(Figure 7A). According to the value of K, we divided the 
emerging simulation results into three regions: Region I, in 
which the cooperative target repression induced by the pair 
of miRNAs is negligible compared with independent target 
repression; Region III, in which the repression of p21 is 
mainly driven by the miRNA cooperative interaction and 
Region II, in which independent and cooperative target 
regulation contributes equally. 

When the strength of cooperativity (K) is weak 
(Figure 7A, Region I), both miR-572 and miR-93 work 
independently to repress p21 gene expression. The 
strongest p21 repression is achieved when both are 
overexpressed (purple triangles). Overexpression of one 
or the other miRNA leads to different levels of target re- 
pression due to their different individual repression 
efficacies. For very high values of K (Figure 7A, Region 
III), p21 repression saturates and enters into an inter- 
dependent mode in which the cooperative interaction over- 
comes the repression by individual miRNAs. p21 silencing 
is achieved when both miRNAs are overexpressed (purple 
triangles). In the intermediate region that describes target 
repression for moderate values of AT (Figure 7A, Region II), 
the performance of the target repression increases grad- 
ually. The most remarkable and steep decrease in p21 ex- 
pression levels can be observed when both miRNAs are 
overexpressed (purple triangles). The results in Region II 
suggest that active modulation of miRNA cooperativity 
can play an important role in determining the efficacy of 
target repression. In line with this, some RNA-binding 
proteins (RBPs) like fragile X mental retardation protein 
(FMRP) or PUmilio-Fem-3-binding factor proteins (PUF) 
can enhance target repression by interacting with the 
miRNA-induced silencing complexes (miRISCs), while 
others (e.g. DND1 and HuR) are reported to counteract 
miRNA-mediated repression (32). We hypothesized that 
the activity of those RBPs in combination with the repres- 
sion by multiple miRNAs here described induces a 
tunable-like target repression. According to this, the mech- 
anism of synergistic target regulation can induce a 
sophisticated regulation of miRNA targets, in which 
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Figure 6. The three regulatory mechanisms proposed, particularized for the miRNA pair miR-93 and miR-572. For each depicted target regulation 
mechanism, the p21 expression level is computed for different combinations of the initial concentrations of miR-93 and miR-572. The black symbols 
stand for four different miRNA initial concentration scenarios: both miRNAs are normally expressed ([1, 1], filled inverted triangle), miR-572 is 
overexpressed ([1, 10], filled diamond), miR-93 is overexpressed ([10, 1], filled circle) and both miRNAs are overexpressed ([10, 10], filled triangle). 
The legend bar represents the p21 expression ranging from basal levels (1) to full repression (0). n.u.: normalized unit. 
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Figure 7. (A) p21 expression levels at different strengths of cooperativity (K) between miR-93 and miR-572. The dashed vertical lines represent 
thresholds of K, for which the miRNA-mediated p21 repression displays different behaviours. Different symbols correspond to different miRNA 
initial concentration scenarios: both miRNAs are normally expressed ([1, 1], filled inverted triangle), miR-572 is overexpressed ([1, 10], filled 
diamonds), miR-93 is overexpressed ([10, 1], filled circles) and both miRNAs are overexpressed ([10, 10], filled triangles). The small plot zooms 
in on the first scenario for illustrating the sigmoid shape of the expression levels of p21 at different K. (B) The role of RBPs in multiple 
miRNA-mediated target repression. Our analysis suggests the activity of RBPs in combination with the repression by multiple miRNAs here 
described can induce a tunable-like target repression, n.u.: normalized unit. 



regulation of RBPs activity shifts the performance of target 
repression from one region to another (Figure 7B). 

The response of p21 to upstream stimulus signals 

We further investigated whether the combination of 
changes in miRNA abundance and the cooperativity 
between miRNAs can affect the features of p21 response 



(p21) to transient stimulation (e.g. upstream signal from 
p53 after response to DNA damage). Subsequently, we 
defined a transient stimulus signal activating the synthesis 
of p21 mRNA (mp21) as a function of [L and r, which 
represent the amplitude and duration of the signal, re- 
spectively (Figure 8 A, left). The upstream signal 
function (S(h,t)) was included into the model. 
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Figure 8. The dynamics of p21 to upstream signals. (A) The sketch of the stimulus signal computed (left) and the p21 response (right). The 
parameters u and x account for the amplitude and duration of the stimulus signal, respectively. (B) The predicted p21 response peak to different 
stimulus signals for different miRNA abundance scenarios. Each plot refers to a set of stimulus signals with the same amplitude range but different 
durations (a 6 [1CF 2 10 1 ] n.u., x = [1, 10, 24] h). Five miRNA abundance scenarios were considered and they are as follows: (i) off (solid red line): no 
expression of any miRNA, (ii) on (black dash-dot line): normal expression of all miRNAs, (iii) on + C (blue circles): normal expression with miRNA 
cooperativity, (iv) on(xlO) (green dotted line): overexpression of all miRNAs and (v) on(xlO) + C (magenta dashed line): overexpression with miRNA 
cooperativity. The grey dashed lines represent the threshold, which is equal to 10% of the basal p21 expression level. (C) The time-series plots of the 
p21 response when the p21 is stimulated by the low (n = 0.1 n.u.) or high (u = 1 n.u.) amplitude and long-lasting (r = lOhr) stimulus signals, n.u.: 
normalized unit. 
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In a series of simulations, we modulated \i and r 
and computed the peak of the p21 response (Figure 8 A 
right). Assuming non-basal synthesis for each miRNA 
{ksyrK-miRKi = i ■ ■ ■ 15) = 0), we simulated the p21 response 
for five different miRNA abundance scenarios by 
modulating their initial concentrations: (i) off. no expres- 
sion of any miRNA (miRi = 1 ■ ■ ■ is(0) = 0, k dass _, niRi j = 



kdeg_cpXij = 0), (ii) on. normal expression of all miRNAs 
(miRi = j... 7J (0) = 1, k class _,„ iRi j = kdeg_cpxij = 0), (iii) 
on + C: normal expression with miRNA cooperativity 

(miRi = 1 ■ ■ ■ /5(0) — 1 , kdass—mUltj = k ass _ m ijn + k ass _ m ijy, 

kdeg_c P xij = k de g_cpxi + k deg _cpxj), (iv) on(xlO): over- 
expression of all miRNAs (miR, = j . . . 15 (0) = 10, 
kdass-jniRij = k c/eg _c P xij = 0) and (v) on(xl0) + C: overex- 
pression with miRNA cooperativity (miR, != / . . . 15 (0) = 

10, k dass m iRij k ass miRi k ass miRp k de g_CpXiJ k de g_ 

CpXi+kaeg_cpxj)- Furthermore, we computed the peaks 
of the p21 response to transient stimulation for differ- 
ent settings of stimulus signals (Figure 8B; see Supple- 
mentary Materials Section 3.2 for details). Moreover, we 
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compared the dynamics of the p21 response for 
long-lasting stimulus signals with low or high amplitude 
(Figure 8C; see Supplementary Materials Section 3.2 for 
details). 

For short stimulation (Figure 8B, x = 1 h), the peak of 
p21 response is positively correlated with the amplitude of 
the stimulus signal, i.e. the stronger the stimulus signal the 
higher the peak. For the complete range of signal ampli- 
tude simulated, the p21 response peak is negatively 
correlated with miRNA abundance. That means the p21 
response peak is highest if miRNAs are not expressed, 
whereas p21 peak is gradually reduced with the increase 
of miRNA abundance (off > on> on(xlO)). In addition, 
miRNA cooperativity further lowers the p21 response 
peak (on> on + C, on(xlO) > on(xlO) +C). However, 
with an increasing duration of the stimulus signal 
(Figure 8B, r = 10 and 24 h), the p21 response peak 
pattern is distorted. For signals with low amplitude 
(H<1), the p21 response peak maintains the same 
pattern (pff> on> on + C> on(xlO) > on(xlO) +C). In 
contrast, for higher amplitude ((i>l) the p21 response 
peak in scenarios with miRNA expression gradually con- 
verges towards a situation of non-miRNA repression. This 
behaviour is due to the consumption of the entire amount 
of available miRNAs (data not shown). Interesting 
enough, higher peaks of p21 response are reached for 
on(xlO) + C than for on(xlO) when the amplitude and 
duration of the stimulus signal are intense and long 
(H>8, r=10h or u>2, r = 24h). This result seems 
contradictory to our expectation of stronger p21 repres- 
sion caused by miRNA cooperation. However, the simu- 
lations of miRNA consumption dynamics indicate that 
under identical miRNA abundance, the model that con- 
siders cooperation shows quicker exhaustion of free 
miRNA, which in turn leads to higher peaks of p21 
response (Supplementary Materials Figure S5A). 

In addition, we defined a hypothetical threshold (10% 
of the basal p21 expression level) above which the p21 
response peak is considered as significant. As shown in 
Figure 8B, the position in which the peak crosses the 
threshold accounts for the minimum signal amplitude 
required to obtain significant p21 response. The increasing 
miRNA abundances shift the cross point towards higher 
signal amplitude required (Figure 8B left). The inclusion 
of miRNA cooperativity further enhances this effect 
(Figure 8B middle and right). On the other hand, the 
increasing duration of the signal reduces the amplitude 
required. These results suggest that higher abundance of 
miRNAs as well as the phenomenon of cooperativity can 
enhance noise buffering by increasing the minimum TF 
activity level required to trigger significant target expres- 
sion (33). 

Moreover, for long stimulus duration increasing abun- 
dance of miRNA and cooperativity not only lowers the 
p21 response peak but also changes the shape of the p21 
temporal pattern (Figure 8C). For low stimulus ampli- 
tude, the p21 response exhibits a rather flat activation, 
with values of p21 in the same interval for many hours, 
while for high signal amplitude the activation of p21 is 
delayed. More simulations with variation in signal 



duration and amplitude are included in Supplementary 
Materials Figure S5B. 

Cooperative miRNA regulation of p21 expression for 
different cellular functions 

We used the model to predict the p21 steady-state levels 
for different cellular functions and to analyse the role of 
miRNA cooperativity. Subsequently, we first extracted the 
list of GO terms associated with the TFs of the miRNAs. 
We then assumed activation of miRNA expression when 
at least one of the miRNA TFs is associated with the 
cellular functions under investigation. Figure 9 (left) 
displays the predicted activity profile of p21 targeting 
miRNAs for different cellular functions. In addition, we 
computed the p21 steady-state levels affected by coopera- 
tive vs non-cooperative miRNA action (Figure 9 right). 
The model predicts high levels of p21 expression during 
DNA damage, DNA repair, senescence and migration, 
while low levels are associated with cell proliferation, 
apoptosis and cell cycle. Moreover, in some of the 
cellular functions p21 steady-state levels are identical for 
both cooperative and non-cooperative miRNA-mediated 
regulation. This can be explained by the fact that the 
miRNAs expressed under these scenarios do not have 
their binding sites in close proximity and therefore do 
not effectively cooperate (Figure 9 right; DNA repair, 
cell migration and DNA damage). In contrast, 
cooperating miRNAs are expressed in other scenarios, 
and thus differences appear in the p21 steady-state levels 
under the cooperative and non-cooperative miRNA 
action (Figure 9 right; senescence and immune response). 
This suggests that cooperating miRNAs can be selectively 
expressed in specific biological scenarios. 

CONCLUSIONS AND DISCUSSION 

Structure of the mathematical model 

So far only a few kinetic models of miRNA target regu- 
lation are proposed. The mathematical model of the 
miRNA-mediated gene silencing developed by Levine 
et al. (34) was continued by Nissan and Parker (35) and 
Zinovyev et al. (36), who developed detailed models 
describing the mechanism of miRNA-induced target re- 
pression. Compared to them, in our model we assumed 
a simplified description for the mechanism of 
miRNA-mediated target repression. However, all these 
models considered only regulation by individual 
miRNAs, while the model presented in this article ad- 
dressed the question of the multiple and cooperative 
miRNAs repressing the same target. Khanin and 
Vinciotti (37) derived a data-driven model using micro- 
array datasets of human mRNAs regulated by miR- 
124 a and used it to quantify the miRNA-mediated effect 
on target mRNA degradation. We followed a similar 
approach using quantitative data to calibrate our model, 
but the model here presented also quantifies the effect of 
miRNA repression in both mRNA and protein levels. Lee 
et al. (38) constructed a network model characterizing 
miR-204 as a tumour suppressor that was used to 
detect/validate 18 gene targets related to tumour 
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Figure 9. Cooperative miRNA regulation of p21 expression in different cellular functions. The statuses (on/off) of miRNAs for specific cellular 
functions are derived from the GO terms of their TFs. The horizontal bars represent the predicted p21 steady-state levels for different cellular 
functions when the miRNA cooperativity is considered (C) or not (no C). The equations and parameter values can be found in Supplementary 
Materials Section 3.3. n.u.: normalized unit. 



progression. In line with this idea, we constructed the 
model based on the p21 regulatory network and used it 
to predict p21 expression for nine different cellular 
functions. 

Regulation of p21 

We here used a novel approach, combining bioinformatics 
algorithms and data from the literature and databases 
with kinetic modelling, to elucidate the mechanisms 
involved in the regulation of the miRNA target hub p21. 
We constructed a map based on TF prediction algorithms, 
published reports and databases, depicting the transcrip- 
tional and post-transcriptional regulation of p21 as well as 
that of its regulating miRNAs. We derived a kinetic model 
to test the hypotheses concerning mechanisms of collective 
miRNA repression and proposed at least three possible 
mechanisms. We computationally detected several 
couples of miRNAs that may cooperate and experimen- 
tally validated the cooperativity between miR-572 and 
miR-93. Our results suggest that these two miRNAs can 
induce cooperative effect to repress the expression of p21 
via the mechanism of synergistic target regulation. 

Through GO analysis, cellular functions were included 
in the map opening the possibility to investigate tissue- 
and development-dependent patterns of p21 regulation. 
Our model predicts high levels of p21 expression during 
DNA damage, DNA repair, senescence and migration, 
while low levels are associated with cell proliferation, 
apoptosis and cell cycle. Our predictions are supported 
by recent experimental publications describing p21 expres- 
sion in the considered cellular functions. Interestingly, en- 
dogenous up-regulation of a number of miRNAs targeting 
p21 can rescue cells from p21-mediated senescence (19). 



These findings are in line with the model predictions that 
the miRNA-mediated p21 repression is poor during 
senescence, but becomes strong when cells re-enter cell 
proliferation. In addition, the model suggests that a low 
level of p21 expression is associated with apoptosis. It is 
known that one of the possible events triggering apoptosis 
initiation is the mitotic catastrophe in which cells fail to 
properly arrest cell cycle (39). Insufficient expression of 
cell cycle arresters like p21 and 14-3-3 a may trigger this 
phenomenon (40,41). Our simulations predict that an 
apoptosis-associated miRNA expression profile can con- 
tribute to the down-regulation of p21 during apoptosis, 
which could facilitate mitotic catastrophe. 

Moreover, our map and kinetic model provide a useful 
means to predict and analyse deregulation in p21 expres- 
sion associated with cancer-related phenotypes. For 
example, we detected a number of feedforward loops 
composed of p21, p21 -targeting miRNAs and their 
common TFs (Supplementary Materials Table SI). The 
existence of feedforward loops has important conse- 
quences for the deregulation of p21 in cancer conditions. 
Interestingly, p21 has dual roles during the cell cycle 
process: under non-stressed conditions, p21 is expressed 
at low levels and promotes cell cycle progression; under 
stress conditions, p21 expression is increased through 
p53-dependent pathways and becomes a cell cycle inhibi- 
tor (42). Our results predict that the tumour suppressor 
p53 could regulate p21 through an incoherent feedforward 
loop, in which p53 directly up-regulates p21 and indirectly 
down-regulates p21 through the p53-promoted transcrip- 
tion of miR-125-5p. In non-small lung cancer cells, 
miR-125a-5p is down-regulated, suggesting that for some 
cancer-associated phenotypes the feedforward loop 
formed by p53, p21 and miR-125a-5p is deactivated 
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(43). Our analysis suggests that the suppression of such 
regulatory loop can favour cancer progression by delaying 
the initiation of p21 -triggered cell cycle arrest after DNA 
damage (Supplementary Materials Figure S7). 

Design principles of networks involving 
miRNA target hubs 

Levine et al. (34) hypothesized that the extent to which an 
miRNA can affect the mRNA levels of its target genes can 
be altered by miRNA-specific parameters, global param- 
eters (controlling mRN A turnover) and the features of the 
cellular function in which repression occurs. Our results 
are in line with their conclusions and indicate that miRNA 
cooperativity or synergy in gene repression may be one of 
the mechanisms by which specificity is gained in the regu- 
lation of different targets by miRNAs. But it also suggests 
that cooperating miRNAs may induce fine-tuning in gene 
expression for cellular functions in which those close 
binding site miRNAs exhibit distinctive regulation. 

Our kinetic model supports the findings of Doench and 
Sharp (6), confirming that repression of certain genes may 
require two or more endogenous miRNAs targeting 
binding sites that are in close proximity. Using a distance 
of 13-35 nt between two seed sites as criterion for allowing 
enhanced repression (7), our computational analysis 
detected many occurrences of binding sites in close prox- 
imity for p21. When integrated and analysed by mathem- 
atical modelling, the regulatory network derived here 
indicates that distinctive regulation, through different 
TFs, associated with different critical cellular functions, 
may be one of the mechanisms by which the cellular 
function determines efficiency of miRNA-mediated repres- 
sion (Figure 10). 

Although we focused on pairwise interactions between 
miRNAs, the prediction of miRNA-binding sites for p21 
uncovered larger putative motifs with three or more 
non-overlapping miRNA-binding sites in close proximity 
(e.g. miR-657, miR-298 and miR-208a; see Figure 4) or 
partially overlapping ones (e.g. miR-654-3p, miR-572 and 
miR-93). These motifs may account for the clustered 
miRNA-binding sites proposed by Rigoutsos et al. (44) 
and suggest the existence of more complex forms of 
nonlinear regression. In many cases, miRNAs can only 
induce mild repression of their targets (4). Based on our 
analysis, the mechanism of interdependent regulation can 
explain the relatively poor ability of overexpression of 
single miRNAs to silence genes, which is observed 
in vitro or in vivo (45). 

Moreover, the mechanism of synergistic target regula- 
tion implies a possible sophisticated regulation of miRNA 
targets. In this mechanism, the strength of miRNA 
cooperativity is determined by the interplay between the 
miRISCs and RBPs (46). 

Our analysis also confirms the enrichment of network 
motifs in the local architecture of p21 regulation, like 
coherent and incoherent feedforward loops containing 
miRNAs (3,47,48). Our map reveals 30 regulatory loops. 
This includes 27 feedforward loops integrated by TFs 
that regulate the expression of miRNAs and p21 
(Supplementary Materials Table SI). This suggests that 



the combination of these regulatory loops with the co- 
operative repression by miRNAs constitutes a novel 
mechanism, by which transient and long-term cell func- 
tion-specific responses are fine-tuned. For example, under 
this assumption the cellular function modulates duration 
and intensity of transient peaks in target protein concen- 
tration predicted for incoherent feedforward loops upon 
step-like TFs activation (49). This regulation can be 
modulated via cooperation with other miRNAs, expressed 
in a cell-function specific manner (Figure 11). 

In our analysis, a kinetic model is used to describe the 
miRNA-mediated repression of targeted genes, although 
we are aware of previous models attempting to elucidate 
the details of the miRNA mechanism of repression 
(35,36). Our results, rather than contradicting their state- 
ments, complement them by explaining an additional level 
of fine-regulation in the miRNA-mediated repression via 
interactions between different miRNAs with common 
targets. Taken together, our results describe three mech- 
anisms, based on the cooperation between miRNAs tar- 
geting the same gene, by which the miRNA target 
specificity could be realized. Our analysis provides a 
feasible explanation for the mild repression induced by 
single miRNAs at basal levels found in vivo (4). The coun- 
terpart of this mild repression would be the ability of 
miRNAs, via cooperation, to fine-tune transient or 
long-term responsiveness and cellular function. 

Approach 

The approach here discussed extends and improves our 
preliminary idea of integrating algorithms for target pre- 
diction and databases to develop more comprehensive 
computational models of miRNA regulatory networks 
(50). Given the existence of up to 834 target hubs in 
Homo sapiens (3), the approach proposed here can be 
used as a template for other miRNA target hub 
networks. Our approach can be used to span a gene regu- 
latory network around a target hub of interest by not only 
implementing miRNA-mediated post-transcriptional 
control mechanisms but also transcriptional regulation 
induced by TFs of the miRNAs and the target hub. The 
procedure described here allows the integration of many 
sources of information describing several levels of regula- 
tion in those networks, including gene activation via TFs, 
miRNA-mediated gene silencing and signalling events. 
This information is organized in detailed and annotated 
regulatory maps. Some resources provide the information 
that is predicted by the use of computational algorithms 
(e.g. miRNA-binding sites prediction). Thus, it may be 
noted that the individual interactions in the p21 regulatory 
map cannot be taken as granted, although we tried to 
implement only those interactions that are likely to 
occur. However, we calculated for each interaction a con- 
fidence score that lies between 0 and 1 and is intuitive for 
experimentalist. Analysis of the confidence scores indi- 
cates that in general terms protein-protein interactions 
and miRNA-target interactions have higher scores 
(respectively, 0.72 and 0.78 on average) because for these 
interactions we only considered highly reliable informa- 
tion, essentially direct interactions that are experimentally 
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Figure 10. Examples of function-specific regulation of p21. The TFs (rectangles) regulate p21 directly and indirectly through miRNAs (parallelo- 
grams): inflammatory response (A) and cell cycle (B). 
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Figure 11. Target hub regulation by cooperative miRNA repression and feedforward loops. The regulatory loops detected in our analysis and 
cooperative miRNA repression can synergize in the regulation of the target. This synergy can generate further fine-tuning in the response. As an 
example, we here display a module which combines an incoherent feedforward loop with transcriptional delay in miRj expression (TF^TgHub, 
TFj^miR!, miRHTgHub) and cooperative target repression ([miR[ AND miRJHTgHub, left). When miR 2 expression is triggered by a different 
transcription factor TF 2 , it can modulate the transient peak and the long-term target expression levels. Our simulations indicate that, beyond a 
threshold in TF 2 -mediated expression of miR 2 , a step-like transcriptional activation of TFi induces a transient peak in target hub expression 
(middle). The intensity of that transient peak is modulated by the levels of TF 2 , which promotes the expression of the second cooperative 
miRNA, miR 2 (right), n.u.: normalized unit. 



validated. We note that there is an immense amount of 
published information on protein-protein interactions, 
while for miRNA-target interactions the purpose of our 
analysis limits us to consider only those whose actual re- 
pressive ability is measured. In line with this, moderately 
high scores were computed for TF-p21 interactions where 
typically at least one publication per interaction is found. 
In contrast, most of the TF-miRNA relationships are 



computationally predicted (investigation of transcrip- 
tional regulation of miRNAs is merely emerging now). 
Therefore, these interactions received rather low scores 
(0.32 average). However, to ensure certain confidence we 
consider only the TF-miRNA interactions predicted by at 
least two independent algorithms. Nevertheless, TF regu- 
lation of miRNAs requires experimental validation in 
most of the cases. Likewise, some of the p21 TFs are 
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putative and demand further investigation. Thus, the 
analysis of our map provides the basis for further experi- 
mental investigations, especially for the non-confirmed 
interactions involved in the regulatory loops. 

The combination of different network motifs in the 
visual map allows the detection of a rich landscape of 
positive and negative interconnected or overlapping regu- 
latory loops. This enables the identification of potential 
network motifs, in particular feedback and feedforward 
loops. The latter ones play a noise buffering role (51). 
Under these conditions, mathematical modelling 
becomes a valuable tool to dissect the behaviour of 
those regulatory loops. 

We show how sections of the regulatory map can be 
translated into ODEs and subsequently be characterized 
with biological information and quantitative data. The 
data should come from perturbation experiments, in 
which expression levels of the miRNAs considered are 
tuned (up- or down-regulated) and the effect on the 
target mRNA and protein levels are quantified. In the 
ideal case, the model should be refined through iterative 
cycles of parameter estimation, identifiability analysis, 
model structure modification and new perturbation ex- 
periments (26,27). 

The deduced regulatory network can subsequently be 
used to investigate the dynamics that target hub expres- 
sion undergoes due to the compound regulation con- 
ducted by dozen(s) of miRNAs and TFs. The proposed 
system of ODEs provides a solid fundament for these in- 
vestigations. Full characterization in terms of model 
parameterization is desired but many times not achievable 
due to the lack of appropriate quantitative data. In this 
case, several analytical methodologies can be used to 
explore the design space of the system (52,53), a simplified 
version of which are used here to investigate the mechan- 
isms of cooperative miRNA-mediated regulation in target 
hubs. 
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